This file analyzes the changes in S/H cell ratio (bleaching and symboint shuffling) in Pocillopora spp.: * 29 colonies were sampled in 2014 (Before ENSO), 2015 (during a 1st peak of heat) and 2016 (after a 2nd peak of heat) * These samples were collected at Uva Island, Gulf of ChiriquÃ, Panama.
# Libraries
library(plyr)
library(tidyverse) # install.packages('tidyverse')
# library(devtools) # install.packages("devtools")
# devtools::install_github("jrcunning/steponeR")
library(steponeR)
library(scales)
library(gridExtra)
library(lme4)
library(lmerTest)
library(multcomp)
library(emmeans)
# Default ggplot settings
ggthe_bw<-theme(plot.background=element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.line = element_line(colour = "black"),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)+
theme_bw()
Get the raw data for all Pocillopora samples applying R.Cunning steponeR function
Pdam.plates <- list.files(path="data", pattern=".csv", full.names=T)
Pdam.plates
## [1] "data/Chi-Pdam_Plate2.csv" "data/Chi-Pdam_Plate3.csv"
## [3] "data/Chi-Pdam_Plate4.csv" "data/Chi-Pdam_Plate5.csv"
## [5] "data/Chi-Pdam_Plate6.csv" "data/Gal-Pdam_Plate13.csv"
## [7] "data/Gal-Pdam_Plate14.csv" "data/Gal-Pdam_Plate15.csv"
## [9] "data/Gal-Pdam_Plate16.csv" "data/Gal1_Dprobes_data.csv"
## [11] "data/Pan-Pdam_Plate10.csv" "data/Pan-Pdam_Plate7.csv"
## [13] "data/Pan-Pdam_Plate8.csv" "data/Pan-Pdam_Plate9.csv"
## [15] "data/Pdam_Plate11-Deep.csv"
Pdam.Out <- steponeR(files=Pdam.plates, target.ratios=c("C.Pdam", "D.Pdam"),
fluor.norm=list(C=3.798, D=0, Pdam=5.416),
copy.number=list(C=9, D=1, Pdam=3),
ploidy=list(C=1, D=1, Pdam=2),
extract=list(C=0.813, D=0.813, Pdam=0.982))
# Target ratio results
Pdam<-Pdam.Out$result
Get colony and year information from sample labels
# Parse sample names, times, treatments, colonies, plates
sample.year <- rbind.fill(lapply(strsplit(as.character(Pdam$Sample.Name), split="_"),
function(X) data.frame(t(X))))
colnames(sample.year) <- c("Spp", "Colony","Year")
Pdam <- cbind(sample.year, Pdam[,-1])
# Keep unique sample ID to check reruns
Pdam$Sample<-paste(Pdam$Year,Pdam$Colony, sep="_")
Pdam$ID<-paste(Pdam$Sample,Pdam$File.Name, sep=".")
Calculate total S/H ratio and D/C ratio and Proportion of Durusdinium in each sample
# 0. If Clade only detected in one technical replicate, set its ratio to NA (becomes zero)
Pdam$C.Pdam[which(Pdam$C.reps==1)] <- NA
Pdam$D.Pdam[which(Pdam$D.reps==1)] <- NA
# 1. Rename cols and make NA=0
colnames(Pdam)[which(colnames(Pdam) %in% c("C.Pdam", "D.Pdam" ))] <- c("C.SH", "D.SH")
Pdam$C.SH[is.na(Pdam$C.SH)] <- 0
Pdam$D.SH[is.na(Pdam$D.SH)] <- 0
Pdam <- Pdam[, c("Colony", "Year", "C.SH", "D.SH","Sample", "ID")]
# 2. Get the ratios and log 10
# Total ratio
Pdam$tot.SH <- Pdam$C.SH + Pdam$D.SH
Pdam$logTot.SH <- log10(Pdam$tot.SH )
# C ratio
Pdam$logC.SH <- log10(Pdam$C.SH)
# D ratio
Pdam$logD.SH <- log10(Pdam$D.SH)
Pdam$logTot.SH[which(Pdam$tot.SH==0)] <- NA
Pdam$logC.SH[which(Pdam$C.SH==0)] <- NA
Pdam$logD.SH[which(Pdam$D.SH==0)] <- NA
# Clade Proportion
# D Proportion
Pdam$D.Prp<-(Pdam$D.SH/Pdam$tot.SH)
# C Proportion
#Pdam$C.Prp<-(Pdam$C.SH/Pdam$tot.SH)
hist(Pdam$D.Prp)
#Dproportion<-Pdam[(Pdam$D.Prp<0.99 & Pdam$D.Prp>0.01),]
#hist(Dproportion$D.Prp)
Note: Most of the colonies are dominated by one symbiont genus. Even proportions of multiple symbionts are rare.
## 1. Check and remove NTC wells
ntc <- Pdam[which(Pdam$Colony=="NTC"), ]
Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(ntc), ])
## 2. Check and remove NoHSratio samples
NoHSratio <- Pdam[which(Pdam$tot.SH==0), ]
Pdam <- droplevels(Pdam[!rownames(Pdam) %in% rownames(NoHSratio), ])
## 3. Chose bw samples ran more than once
ReRunA <- Pdam[duplicated(Pdam$Sample),]
n_RunA <- data.frame(table(Pdam$Sample))
colnames(n_RunA)<-c("Sample","RanA")
Pdam<-join(Pdam, n_RunA, type = "left")
DuplicatesA <- Pdam[(Pdam$RanA>1),]
# Remove bad replicates
ToRem1<-read.csv("ToRemove10-28-16.csv") # 11/1/16
Pdam<-Pdam[!(Pdam$ID %in% ToRem1$ID),]
n_RunB <- data.frame(table(Pdam$Sample))
ReRunB <- Pdam[duplicated(Pdam$Sample),]
colnames(n_RunB)<-c("Sample","RanB")
Pdam<-join(Pdam, n_RunB, type = "left")
# List of dupplicated samples, should have 0 rows now
DuplicatesB <- Pdam[(Pdam$RanB>1),]
# 4. Check and remove samples with high SD or late amplification
StDe1.5 <- Pdam[which(Pdam$Pdam.CT.sd>1.5), ]
Pdam<-Pdam[!(Pdam$ID %in% StDe1.5$ID),]
# End of qPCR data cleaning
Get collection information and select only samples from Uva Island that were sampled in 2014, 2015 and 2016
# 1. Labeling and Factors
## Import treatment info from Field file and mortality / ITS2 / ORF data
Info<-read.csv("SampleInfo.csv", header=TRUE)
Mortality<-read.csv("MortalityData.csv", header=TRUE)
Mortality$ITS_2<-factor(Mortality$ITS_2, level=c("C1d", "C1b-c", "D1"))
## Merge SH ratios and sample info
Pdam<-join(Pdam, Info, by = "Sample",
type = "left", match = "all")
#summary(Pdam)
# 2. Uva Data
# Write how many years a colony was sampled
Years <- data.frame(table(Pdam$Colony))
colnames(Years)<-c("Colony","Times")
Pdam<-join(Pdam, Years, type = "left")
# Select data only from Uva Island
data.Uva<-Pdam[which(Pdam$Location=="Uva Island"), ]
# Select colonies that were sampled three times (2014, 2015 and 2016)
data.Uva<-data.Uva[which(data.Uva$Times>2),]
# Remove "Ross" samples (colony in deep location (~50feet), not Uva reef)
# and 259 as well as 264 that were mislabeled in at least one of the sampling points
data.Uva<-filter(data.Uva, Colony!="Ross")
data.Uva<-filter(data.Uva, Colony!="259-1")
data.Uva<-filter(data.Uva, Colony!="264-1")
data.Uva<-data.Uva[,c("Colony","Year.2","C.SH","D.SH","tot.SH",
"logTot.SH","D.Prp", "logC.SH", "logD.SH",
"Year","Sample", "Reef", "Spp", "Depht")]
data.Uva$Year<-factor(data.Uva$Year, levels=c("14", "15", "16"))
data.Uva$Colony<-factor(data.Uva$Colony)
#summary(data.Uva)
#3. Merge SH ratios and moratlity info
data.Uva<-join(data.Uva, Mortality, by = "Sample",
type = "left", match = "all")
summary(data.Uva)
## Colony Year.2 C.SH D.SH
## 256-1 : 3 Min. :2014 Min. :0.000000 Min. :0.00000
## 257-1 : 3 1st Qu.:2014 1st Qu.:0.000000 1st Qu.:0.00000
## 260-1 : 3 Median :2015 Median :0.003005 Median :0.01442
## 262-1 : 3 Mean :2015 Mean :0.019174 Mean :0.02071
## 278-1 : 3 3rd Qu.:2016 3rd Qu.:0.011841 3rd Qu.:0.03001
## 284-1 : 3 Max. :2016 Max. :0.690021 Max. :0.18875
## (Other):69
## tot.SH logTot.SH D.Prp logC.SH
## Min. :0.00146 Min. :-2.8356 Min. :0.0000 Min. :-5.0627
## 1st Qu.:0.01294 1st Qu.:-1.8888 1st Qu.:0.0000 1st Qu.:-2.5001
## Median :0.02266 Median :-1.6448 Median :0.8620 Median :-2.0971
## Mean :0.03988 Mean :-1.6539 Mean :0.5411 Mean :-2.1808
## 3rd Qu.:0.03579 3rd Qu.:-1.4463 3rd Qu.:1.0000 3rd Qu.:-1.7015
## Max. :0.69002 Max. :-0.1611 Max. :1.0000 Max. :-0.1611
## NA's :32
## logD.SH Year Sample Reef
## Min. :-5.7262 14:29 Length:87 Coral Collection Point:36
## 1st Qu.:-1.8381 15:29 Class :character Shallow Millepora :18
## Median :-1.6245 16:29 Mode :character M.intrincata Shallow :12
## Mean :-2.0059 Ridley Rock Shallow : 9
## 3rd Qu.:-1.4640 Gardineroseris city : 6
## Max. :-0.7241 Bw Gard City and 4x5 : 3
## NA's :29 (Other) : 3
## Spp Depht Mortality ID_Glynn ID_Palacio
## P.dam:51 Min. : 9.00 Almost total: 7 P. cap:12 P. dam:24
## P.eff: 0 1st Qu.:15.00 No :71 P. dam:27 P. eyd: 6
## P.ele:36 Median :17.00 Partial : 9 P. ver:48 P. ver:57
## P.eyd: 0 Mean :17.79
## P.mea: 0 3rd Qu.:20.00
## P.ver: 0 Max. :32.00
##
## ORF ITS_2
## Type_1:69 C1d :18
## Type_3:18 C1b-c:20
## D1 :49
##
##
##
##
data.Uva$Mortality<-factor(data.Uva$Mortality, levels =
c("No", "Partial", "Almost total"))
# Summary ID *Pocillopora* by Glynn
Summary_ID_Glynn<-data.Uva %>%
group_by(Year.2, ORF, ID_Glynn) %>%
dplyr::summarise(total.count=n())
Summary_ID_Glynn
# Summary ID *Pocillopora* by Palacio
Summary_ID_Glynn<-data.Uva %>%
group_by(Year.2, ORF, ID_Palacio) %>%
dplyr::summarise(total.count=n())
Summary_ID_Glynn
# Add a column for the genera detected in each sample using qPCR
data.Uva$Community<-NULL
data.Uva$Community[which(data.Uva$D.Prp>=0)] <- "Mixed"
data.Uva$Community[which(data.Uva$D.Prp==0)] <- "Cladocopium-only"
data.Uva$Community[which(data.Uva$D.Prp==1)] <- "Durusdinium-only"
data.Uva$Community<-factor(as.character(data.Uva$Community),
levels=c("Cladocopium-only", "Mixed", "Durusdinium-only"))
# Summary genera *Pocillopora*
Summary_Community<-data.Uva %>%
group_by(Year.2, Community) %>%
dplyr::summarise(total.count=n())
Summary_Community
# Summary genera detected by ORF type
Summary_Community<-data.Uva %>%
group_by(ORF, Year.2, Community) %>%
dplyr::summarise(total.count=n())
Summary_Community
In general, Pocillopora colonies are dominated by one Symbiodiniaceae type that accounts for more than 80-90 of the symbiont community. However, sometimes intermediate proportions are common when the symbiont community is transitioning from being dominated by one algal type to other (e.g. during heat stress).
For this data set, 2 of the colonies sampled in 2014 had relativelly even proportions of Cladocopium_C1bc and Durusdinium glynnii (Colony 267_2014, D proportion = 0.47; Colony 269_2014, D proportion = 0.50). This made difficult their classification into one domminant type for statistical testing and further plotting. Since these 2 colonies were later dominated by Durusdinim glynnii in 2015, we decided to include them in the D-dominated classification, eventhoug in 2014 there were not such a dominance (See figure 2a).
# Add a column for the dominant genus
data.Uva$DominantCommunity<-NULL
data.Uva$DominantCommunity[which(data.Uva$D.Prp>=0.45)] <- "Durusdinium-dominated"
data.Uva$DominantCommunity[which(data.Uva$D.Prp<0.45)] <- "Cladocopium-dominated"
data.Uva$DominantCommunity<-factor(as.character(data.Uva$DominantCommunity),
levels=c("Cladocopium-dominated","Durusdinium-dominated"))
# Summary dominant genera *Pocillopora*
Summary_Community2<-data.Uva %>%
group_by(Year.2, DominantCommunity) %>%
dplyr::summarise(total.count=n())
Summary_Community2
# Summary dominant genera by ORF type
Summary_Community3<-data.Uva %>%
group_by(ORF, Year.2, DominantCommunity) %>%
dplyr::summarise(total.count=n())
Summary_Community3
Set initial (2014) and final (2016) communities for each colony
data.Uva$ITS_2<-factor(data.Uva$ITS_2, level=c("C1b-c", "C1d", "D1"))
# Subset Aug 2014, Aug 2015 and April 2016 samples
D2014.data<-subset(data.Uva, Year.2==2014)
D2014.data<-D2014.data[,c("Colony","D.Prp","Community","DominantCommunity", "ITS_2")]
colnames(D2014.data)<-(c("Colony","D14.D.Prp","D14Community", "D14Dominant", "ITS2_14"))
data.Uva<-left_join(data.Uva, D2014.data, by=c("Colony")) # Add initial community to all samples
D2016.data<-subset(data.Uva, Year.2==2016) # D1a under control temperature
D2016.data<-D2016.data[, c("Colony","D.Prp","Community","DominantCommunity", "ITS_2")]
colnames(D2016.data)<-(c("Colony","D16.D.Prp","D16Community", "D16Dominant", "ITS2_16"))
data.Uva<-left_join(data.Uva, D2016.data, by=c("Colony")) # Add final community to all samples
To test the statistical significance of changes in the number of colonies dominated by Durusdinium versus Cladocopium we considered inclussion and exclussion of these 2 initially even colonies, and to place them in the D-dominated or C-dominated categories in 2014.
Input =("
Year C D
2014 15 14
2016 8 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.1064
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.8336155 9.7658946
## sample estimates:
## odds ratio
## 2.760908
Input =("
Year C D
2014 16 13
2016 8 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.06099
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.9556883 11.2396870
## sample estimates:
## odds ratio
## 3.162767
Input =("
Year C D
2014 15 12
2016 8 19
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.09778
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.848568 10.659493
## sample estimates:
## odds ratio
## 2.906792
Input =("
Year C D
2014 17 12
2016 8 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.03295
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.094473 13.006776
## sample estimates:
## odds ratio
## 3.629777
Histo_Dprop <- ggplot(data.Uva, aes (D.Prp , fill=(D14Community))) + ggthe_bw +
geom_histogram(binwidth = 0.05, boundary=0) +
facet_grid(~Year.2, labeller = labeller(Year.2=c(
`2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~(D/H)/(S/H)),
breaks=(seq(0,1,by=0.1))) +
scale_y_continuous((name= "Number of colonies"),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
Histo_Dprop
Traject_Proportion_Facet <- ggplot (data.Uva) +
ggthe_bw +
geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+
facet_wrap(~D14Community, labeller = labeller(
D14Community=c(`Cladocopium-only` = "Cladocopium-only (2014)",
`Mixed` = "Mixed (2014)",
`Durusdinium-only` = "Durusdinium-only (2014)"))) +
geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3,
height = 0.05, width = 0.15) +
#labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
theme(legend.position="none")
Traject_Proportion_Facet
Histo2014_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(D14Community))) +
ggthe_bw + coord_flip()+ ggtitle("a.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~((D/H)/(S/H))),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2014_Dprop
Histo2016_Dprop <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(D14Community))) +
ggthe_bw + coord_flip()+ ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2016_Dprop
Traject_Proportion<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp , colour=(D14Community), group=(Colony))) + ggthe_bw + ggtitle("b.") +
geom_line(linetype=2) +
geom_jitter(height = 0.01, width = 0.25, alpha=0.4) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial symbiont community (August 2014)",
labels=c("Cladocopium-only", "Mixed", "Durusdinium-only")) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.01,0.01)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Traject_Proportion
Figure_2<-grid.arrange(Histo2014_Dprop,Traject_Proportion, Histo2016_Dprop,
nrow=1, widths=c(3/9, 3/9, 3/9))
# ggsave(file="Outputs/Figure_2qPCR.svg", plot=Figure_2, width=7, height=4)
SH_2014Co<-lmerTest::lmer(logTot.SH ~ Year * D14Community + (1|Colony), data=data.Uva)
summary(SH_2014Co)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 83.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6538 -0.5202 -0.1347 0.3700 2.9430
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.03036 0.1742
## Residual 0.10764 0.3281
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.6402 0.1238 71.1153 -13.246
## Year15 -0.4132 0.1547 52.0000 -2.672
## Year16 -0.5536 0.1547 52.0000 -3.579
## D14CommunityMixed -0.0137 0.1638 71.1153 -0.084
## D14CommunityDurusdinium-only 0.1232 0.1805 71.1153 0.682
## Year15:D14CommunityMixed 0.7731 0.2046 52.0000 3.778
## Year16:D14CommunityMixed 0.6960 0.2046 52.0000 3.402
## Year15:D14CommunityDurusdinium-only 0.4584 0.2255 52.0000 2.033
## Year16:D14CommunityDurusdinium-only 0.3869 0.2255 52.0000 1.716
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.010053 *
## Year16 0.000756 ***
## D14CommunityMixed 0.933599
## D14CommunityDurusdinium-only 0.497283
## Year15:D14CommunityMixed 0.000408 ***
## Year16:D14CommunityMixed 0.001294 **
## Year15:D14CommunityDurusdinium-only 0.047152 *
## Year16:D14CommunityDurusdinium-only 0.092125 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 D14CmM D14CD- Y15:D14CM Y16:D14CM Y15:D14CD
## Year15 -0.624
## Year16 -0.624 0.500
## D14CmmntyMx -0.756 0.472 0.472
## D14CmmntyD- -0.686 0.428 0.428 0.519
## Yr15:D14CmM 0.472 -0.756 -0.378 -0.624 -0.324
## Yr16:D14CmM 0.472 -0.378 -0.756 -0.624 -0.324 0.500
## Yr15:D14CD- 0.428 -0.686 -0.343 -0.324 -0.624 0.519 0.259
## Yr16:D14CD- 0.428 -0.343 -0.686 -0.324 -0.624 0.259 0.519 0.500
#step(SH_2014Co)
SH.emmc<-emmeans(SH_2014Co, ~Year*D14Community)
SH_groups<-multcomp::cld(SH.emmc)
SH_groups<-SH_groups[order(
SH_groups$Year, SH_groups$D14Community),]
SH_groups
plot(emmeans(SH_2014Co, ~Year|~D14Community), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~D14Community)
SH_2014Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(D14Community))) + ggthe_bw +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position=position_dodge(width=0.95))+
stat_summary(fun.data = "mean_cl_boot",geom = "point",
position=position_dodge(width=0.95))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Symbionts detectected in 2014",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.5,0.1),
strip.background =element_rect(fill="white")) +
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_discrete(name="Year")
SH_2014Com
SH_Community<-lmerTest::lmer(logTot.SH ~ Year * Community + (1|Colony), data=data.Uva)
summary(SH_Community)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * Community + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 91.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3973 -0.5171 -0.0987 0.3905 4.0587
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.02404 0.1550
## Residual 0.12526 0.3539
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.60009 0.12804 76.71004 -12.497 < 2e-16
## Year15 -0.20778 0.15540 49.73821 -1.337 0.187280
## Year16 -0.71036 0.17979 47.25498 -3.951 0.000258
## CommunityMixed -0.09453 0.16829 77.96039 -0.562 0.575910
## CommunityDurusdinium-only 0.09882 0.18645 77.17622 0.530 0.597641
## Year15:CommunityMixed 0.58170 0.24003 54.72489 2.423 0.018708
## Year16:CommunityMixed 0.81742 0.24650 54.96060 3.316 0.001622
## Year15:CommunityDurusdinium-only 0.27313 0.22958 47.93526 1.190 0.240035
## Year16:CommunityDurusdinium-only 0.60594 0.24059 49.84010 2.519 0.015039
##
## (Intercept) ***
## Year15
## Year16 ***
## CommunityMixed
## CommunityDurusdinium-only
## Year15:CommunityMixed *
## Year16:CommunityMixed **
## Year15:CommunityDurusdinium-only
## Year16:CommunityDurusdinium-only *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmntM CmmnD- Y15:CM Y16:CM Y15:CD
## Year15 -0.733
## Year16 -0.602 0.497
## CommuntyMxd -0.754 0.576 0.457
## CmmntyDrsd- -0.686 0.504 0.414 0.524
## Yr15:CmmntM 0.472 -0.658 -0.322 -0.628 -0.326
## Yr16:CmmntM 0.459 -0.365 -0.740 -0.608 -0.318 0.412
## Yr15:CmmnD- 0.496 -0.677 -0.337 -0.383 -0.717 0.440 0.243
## Yr16:CmmnD- 0.451 -0.369 -0.747 -0.326 -0.682 0.240 0.543 0.541
#step(SH_Community)
SH.emmc<-emmeans(SH_Community, ~Year*Community)
SH_groups<-multcomp::cld(SH.emmc)
SH_groups<-SH_groups[order(
SH_groups$Year, SH_groups$Community),]
SH_groups
#write.csv(SH_groups,"SH_groupsB.csv")
plot(emmeans(SH_Community, ~Year|~Community), comparisons = TRUE) +
coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~Community)
SH_Com <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position=position_dodge(width=0.95))+
stat_summary(fun.data = "mean_cl_boot",geom = "point",
position=position_dodge(width=0.95))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.95)) +
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Symbionts detectec in 2014",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.5,0.1),
strip.background =element_rect(fill="white")) +
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_discrete(name="Year")
SH_Com
anova(SH_2014Co, SH_Community)
SH_Communities<-lmerTest::lmer(logTot.SH ~ Year * D14Community * Community + (1|Colony), data=data.Uva)
summary(SH_Communities)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * D14Community * Community + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 81.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.73407 -0.48463 -0.06579 0.35339 2.91935
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.02179 0.1476
## Residual 0.11026 0.3320
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.6402 0.1211 70.1550 -13.541
## Year15 -0.4132 0.1565 47.2363 -2.640
## Year16 -0.7093 0.1688 50.6884 -4.202
## D14CommunityMixed -0.7146 0.3265 73.0936 -2.189
## D14CommunityDurusdinium-only -0.6296 0.3928 72.9975 -1.603
## CommunityMixed 0.7009 0.2845 70.5547 2.464
## CommunityDurusdinium-only 0.7527 0.3508 70.9827 2.146
## Year15:D14CommunityMixed 1.5019 0.3806 64.8657 3.946
## Year16:D14CommunityMixed 0.8259 0.2396 53.9251 3.447
## Year15:D14CommunityDurusdinium-only 1.2674 0.5110 61.4046 2.480
## Year16:D14CommunityDurusdinium-only 0.5426 0.2368 48.9752 2.292
## Year15:CommunityMixed -0.7379 0.3653 70.6556 -2.020
## Year15:CommunityDurusdinium-only -0.8090 0.4572 64.7392 -1.770
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.011206 *
## Year16 0.000107 ***
## D14CommunityMixed 0.031805 *
## D14CommunityDurusdinium-only 0.113266
## CommunityMixed 0.016178 *
## CommunityDurusdinium-only 0.035332 *
## Year15:D14CommunityMixed 0.000198 ***
## Year16:D14CommunityMixed 0.001108 **
## Year15:D14CommunityDurusdinium-only 0.015876 *
## Year16:D14CommunityDurusdinium-only 0.026258 *
## Year15:CommunityMixed 0.047162 *
## Year15:CommunityDurusdinium-only 0.081500 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 14 columns / coefficients
#step(SH_Communities)
SH.emmc<-emmeans(SH_Communities, ~Year*D14Community*Community)
SH_groups<-cld(SH.emmc)
SH_groups<-as.data.frame(SH_groups[complete.cases(SH_groups),])
SH_groups<-SH_groups[order(SH_groups$Year,
SH_groups$D14Community,
SH_groups$Community),]
SH_groups
#write.csv(SH_groups,"SH_Communities_multicomp.csv")
SH_Comties <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(Community))) + ggthe_bw+
facet_grid(~D14Community, labeller = labeller(D14Community=c(
`Cladocopium-only` = "Cladocopium-only (2014)",
`Mixed` = "Mixed (2014)",
`Durusdinium-only` = "Durusdinium-only (2014)"))) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.2),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-3, -0.8),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(14, 15, 16),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
SH_Comties
#ggsave(file="Figure_3.svg", plot=SH_Comties, width=6, height=4)
CH_Com <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(DominantCommunity))) + ggthe_bw+ ggtitle("a.")+
#geom_point()+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Cladocopium", "Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.8),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-3, -0.7),
name="\n Cladocopium to host cell ratio (C/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" ", " ", " "))
#CH_Com
CH_Com2<-CH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
`Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
`Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))
CH<-lmer(logC.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
summary(CH)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * Community * D14Dominant + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 113.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.70311 -0.36865 -0.09556 0.50582 1.89250
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.1182 0.3438
## Residual 0.4032 0.6350
## Number of obs: 55, groups: Colony, 21
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -1.54068 0.23810 44.78323 -6.471
## Year15 -0.25625 0.28031 26.58185 -0.914
## Year16 -0.73614 0.32339 24.71137 -2.276
## CommunityMixed -0.32887 0.37022 44.99572 -0.888
## D14DominantDurusdinium-dominated -1.09559 0.41353 43.96389 -2.649
## Year15:CommunityMixed 0.59388 0.61658 32.70059 0.963
## Year16:CommunityMixed -0.03386 0.50241 28.81910 -0.067
## Year15:D14DominantDurusdinium-dominated -0.54917 0.68606 30.91184 -0.800
## Year16:D14DominantDurusdinium-dominated 0.88798 0.66184 31.43058 1.342
## Pr(>|t|)
## (Intercept) 6.36e-08 ***
## Year15 0.3688
## Year16 0.0318 *
## CommunityMixed 0.3791
## D14DominantDurusdinium-dominated 0.0112 *
## Year15:CommunityMixed 0.3425
## Year16:CommunityMixed 0.9467
## Year15:D14DominantDurusdinium-dominated 0.4296
## Year16:D14DominantDurusdinium-dominated 0.1893
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmntM D14DD- Y15:CM Y16:CM Y15:D1
## Year15 -0.715
## Year16 -0.579 0.492
## CommuntyMxd -0.622 0.502 0.369
## D14DmnntDr- -0.019 -0.038 0.003 -0.537
## Yr15:CmmntM 0.318 -0.474 -0.224 -0.512 0.275
## Yr16:CmmntM 0.400 -0.331 -0.659 -0.643 0.345 0.357
## Yr15:D14DD- 0.006 0.017 0.000 0.255 -0.468 -0.705 -0.185
## Yr16:D14DD- -0.021 0.011 0.011 0.308 -0.509 -0.162 -0.437 0.306
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
#step(CH)
CH.emmc<-emmeans(CH, ~Year*Community* D14Dominant)
CH_groups<-cld(CH.emmc)
CH_groups<-as.data.frame(CH_groups[complete.cases(CH_groups),])
CH_groups<-CH_groups[order(CH_groups$Year,CH_groups$D14Dominant, CH_groups$Community),]
CH_groups
#write.csv(CH_groups,"CH_groups.csv")
DH_Com <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(Community))) +
ggthe_bw+ ggtitle("b.")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3), alpha=0.5)+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Detected symbionts",
labels=c("Mixed", "Durusdinium"))+
theme(legend.position=c(0.82,0.2),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-6, -0.7),
name="\n Durusdinium to host cell ratio (D/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
# DH_Com
DH_Com2<-DH_Com + facet_grid(~D14Dominant, labeller = labeller(D14Dominant=c(
`Cladocopium-dominated` = "Cladocopium- \n dominated (2014)",
`Durusdinium-dominated` = "Durusdinium- \n dominated (2014)")))
DH_LM<-lmer(logD.SH ~ Year* Community* D14Dominant + (1|Colony), data=data.Uva)
summary(DH_LM)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * Community * D14Dominant + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 66.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0953 -0.5756 0.0079 0.3844 3.6310
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.0000 0.0000
## Residual 0.1687 0.4107
## Number of obs: 58, groups: Colony, 22
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) -4.86573 0.16766
## Year15 0.07049 0.33532
## Year16 3.21459 0.23711
## CommunityDurusdinium-only 0.27806 0.50992
## D14DominantDurusdinium-dominated 3.15614 0.23711
## Year15:CommunityDurusdinium-only -0.35506 0.32897
## Year16:CommunityDurusdinium-only -0.22094 0.38416
## Year15:D14DominantDurusdinium-dominated 0.36688 0.42745
## Year16:D14DominantDurusdinium-dominated -3.09697 0.41068
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated -0.08557 0.45916
## df t value
## (Intercept) 48.00000 -29.021
## Year15 48.00000 0.210
## Year16 48.00000 13.557
## CommunityDurusdinium-only 48.00000 0.545
## D14DominantDurusdinium-dominated 48.00000 13.311
## Year15:CommunityDurusdinium-only 48.00000 -1.079
## Year16:CommunityDurusdinium-only 48.00000 -0.575
## Year15:D14DominantDurusdinium-dominated 48.00000 0.858
## Year16:D14DominantDurusdinium-dominated 48.00000 -7.541
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated 48.00000 -0.186
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year15 0.834
## Year16 < 2e-16 ***
## CommunityDurusdinium-only 0.588
## D14DominantDurusdinium-dominated < 2e-16 ***
## Year15:CommunityDurusdinium-only 0.286
## Year16:CommunityDurusdinium-only 0.568
## Year15:D14DominantDurusdinium-dominated 0.395
## Year16:D14DominantDurusdinium-dominated 1.09e-09 ***
## CommunityDurusdinium-only:D14DominantDurusdinium-dominated 0.853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 CmmnD- D14DD- Y15:CD Y16:CD Y15:D1 Y16:D1
## Year15 -0.500
## Year16 -0.707 0.354
## CmmntyDrsd- 0.000 0.000 -0.232
## D14DmnntDr- -0.707 0.354 0.500 -0.232
## Yr15:CmmnD- 0.000 0.000 0.000 -0.293 0.360
## Yr16:CmmnD- 0.000 0.000 0.000 -0.753 0.309 0.389
## Yr15:D14DD- 0.392 -0.784 -0.277 0.129 -0.555 -0.500 -0.171
## Yr16:D14DD- 0.408 -0.204 -0.577 0.671 -0.577 -0.208 -0.713 0.320
## CmD-:D14DD- 0.000 0.000 0.258 -0.900 0.000 0.000 0.558 0.000 -0.596
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#step(DH_LM)
DH.emmc<-emmeans(DH_LM, ~Year*Community* D14Dominant)
DH_groups<-cld(DH.emmc)
DH_groups<-as.data.frame(DH_groups[complete.cases(DH_groups),])
DH_groups<-DH_groups[order(DH_groups$Year,DH_groups$D14Dominant, DH_groups$Community),]
DH_groups
#write.csv(DH_groups,"DH_groups.csv")
Figure_4<-grid.arrange(CH_Com2, DH_Com2,nrow=2)
# ggsave(file="Outputs/Figure_4.svg", plot=Figure_4, width=7, height=4)
#data.Uva$Counts<-paste(data.Uva$DominantCommunity, data.Uva$ORFCommunity, sep = "_")
data.Uva$ORFCommunity<-paste(data.Uva$ORF, data.Uva$Community, sep = "_")
data.Uva$ORFCommunity<-factor(data.Uva$ORFCommunity,
levels = c("Type_3_Cladocopium-only",
"Type_1_Cladocopium-only",
"Type_1_Mixed",
"Type-1_Durusdinium-only"))
data.Uva$ORFCommunity2<-paste(data.Uva$ORF,
data.Uva$DominantCommunity, sep = "_")
data.Uva$ORFCommunity2<-as.factor(data.Uva$ORFCommunity2)
data.Uva$DominantITS2[which(
(data.Uva$DominantCommunity=="Durusdinium-dominated") &
(data.Uva$ORF=="Type_1")
)] <- "D1-dominated"
data.Uva$DominantITS2[which(
(data.Uva$DominantCommunity=="Cladocopium-dominated") &
(data.Uva$ORF=="Type_1")
)] <- "C1bc-dominated"
data.Uva$DominantITS2[which(
(data.Uva$DominantCommunity=="Cladocopium-dominated") &
(data.Uva$ORF=="Type_3")
)] <- "C1d-dominated"
data.Uva$DominantITS2<-as.factor(data.Uva$DominantITS2)
Summary_Community4<-data.Uva %>%
group_by(Year.2, DominantITS2) %>%
dplyr::summarise(total.count=n())
Summary_Community4
Summary_Community5<-data.Uva %>%
group_by(ORF, Year.2, DominantITS2) %>%
dplyr::summarise(total.count=n())
Summary_Community5
# Only type 1 colonies (n=23), p-value = 0.03514
Input =("
Year C D
2014 9 14
2016 2 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.03514
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.105302 70.563784
## sample estimates:
## odds ratio
## 6.4764
# Type 1 and type 3 colonies: p-value = 0.1064
Input =("
Year C D
2014 15 14
2016 8 21
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.1064
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.8336155 9.7658946
## sample estimates:
## odds ratio
## 2.760908
data.Uva$ITS2_14<-factor(data.Uva$ITS2_14, levels = c("C1d", "C1b-c", "D1"))
data.Uva$ITS2_16<-factor(data.Uva$ITS2_16, levels = c("C1d", "C1b-c", "D1"))
Histo_Dprop_ITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(ITS2_14)), shape=(D14Community)) + ggthe_bw +
geom_histogram(binwidth = 0.05, boundary=0) +
facet_grid(~Year.2, labeller = labeller(Year.2=c(
`2014` = "August 2014", `2015` = "August 2015", `2016` = "April 2016"))) +
scale_fill_manual(values=c("blue3", "purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial ITS2 community (August 2014)",
labels=c("Type 3 + C1d", "Type 1 + C1b-c", "Type 1 + D1")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~(D/H)/(S/H)),
breaks=(seq(0,1,by=0.1))) +
scale_y_continuous((name= "Number of colonies"),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
Histo_Dprop_ITS2 +coord_flip()
Histo_Dprop_ORF <- ggplot(data.Uva, aes (D.Prp , fill=(ORF)), shape=(D14Community)) + ggthe_bw +
geom_histogram(binwidth = 0.05, boundary=0) +
facet_grid(~Year.2, labeller = labeller(Year.2=c(
`2014` = "August 2014",
`2015` = "August 2015", `2016` = "April 2016"))) +
scale_fill_manual(values=c("black","gray"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "ORF type",
labels=c("Type 3", "Type 1")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~(D/H)/(S/H)),
breaks=(seq(0,1,by=0.1))) +
scale_y_continuous((name= "Number of colonies"),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
Histo_Dprop_ORF +coord_flip()
# By ITS2 type
Traject_Proportion_FacetITS2 <- ggplot (data.Uva) +
ggthe_bw +
geom_line(aes (x=Year.2, D.Prp, colour=factor(Colony)))+
facet_wrap(~ITS2_14) +
geom_jitter(aes(Year.2, D.Prp, colour=factor(Colony)), alpha=0.3,
height = 0.05, width = 0.15) +
#labs(x = "Año del muestreo", y = "Durusdinium Proportion") +
theme(legend.position="none")
Traject_Proportion_FacetITS2
# By Depth
Traject_Proportion_FacetORF <- ggplot (data.Uva) + ggthe_bw +
geom_line(aes (x=Year.2, D.Prp, group=Colony, linetype=(ORF)))+
geom_jitter(aes(Year.2, D.Prp, colour=(Depht), shape=ORF),
height = 0.00, width = 0.25, size=3, alpha=0.7) +
scale_colour_gradient(low = "red", high = "blue")
Traject_Proportion_FacetORF
Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(ITS2_14))) +
ggthe_bw + ggtitle("a.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Dominant ITS2 community (August 2014)",
labels=c("C1d", "C1b-c", "D.glynnii")) +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2014_DpropITS2
Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(ITS2_16))) +
ggthe_bw + ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Dominant ITS2 community (April 2016)",
labels=c("C1d", "C1b-c", "D.glynnii")) +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2016_DpropITS2
Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp ,
fill=ITS2_14, shape=Mortality, group=(Colony)))+
ggthe_bw + ggtitle("b.") +
geom_line(aes(colour=ITS2_14), linetype=2) +
geom_jitter(height = 0, width = 0.25, alpha=0.4) +
scale_colour_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial ITS2 community (August 2014)",
labels=c("C1d", "C1b-c", "D.glynnii")) +
scale_fill_manual(values=c("blue3","purple", "red3"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Initial ITS2 community (August 2014)",
labels=c("C1d", "C1b-c", "D.glynnii")) +
scale_shape_manual(values=c(21, 22, 24),
labels=c("No", "Partial", "Almost total"))+
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.01,0.01)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Traject_ProportionITS2
Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(2/9, 5/9, 2/9))
# ggsave(file="Outputs/Figure_ITS2.svg", plot=Figure_2_ITS2, width=6.5, height=5)
Histo2014_DpropITS2b<- Histo2014_DpropITS2 + facet_grid(ORF~.)
Histo2016_DpropITS2b<- Histo2016_DpropITS2 + facet_grid(ORF~.)
Traject_ProportionITS2b<- Traject_ProportionITS2 + facet_grid(ORF~.)
Figure_2_ITS2b<-grid.arrange(Histo2014_DpropITS2b, Traject_ProportionITS2b,Histo2016_DpropITS2b,
nrow=1, widths=c(3/9, 3/9, 3/9))
# ggsave(file="Outputs/Figure_2_ITS2b.svg", plot=Figure_2_ITS2b, width=7, height=4)
Histo2014_DpropORF <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(ORFCommunity))) +
ggthe_bw + coord_flip()+ ggtitle("a.")+
geom_histogram(binwidth = 0.05, boundary=0, colour="black") +
scale_fill_manual(values=c("gray", "white","snow", "mintcream"),
guide = guide_legend(title.position = "top", nrow = 2),
name = "")+
#labels=c("Type 3", "Type 1")) +
scale_x_continuous(#name= expression(
#italic(Durusdinium)~proportion~((D/H)/(S/H))),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
limits = c(0,16),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))
#Histo2014_DpropORF
Histo2016_DpropORF <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(ORFCommunity))) +
ggthe_bw + coord_flip()+ ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0,
colour="black") +
scale_fill_manual(values=c("gray", "white","snow", "mintcream"),
guide = guide_legend(title.position = "top", nrow = 2),
name = "")+
#labels=c("Type 3", "Type 1")) +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,16),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA),
panel.border = element_blank(),
axis.line = element_line(colour = "black"))
#Histo2016_DpropORF
Traject_ProportionORF_a<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp))+
ggthe_bw + ggtitle("b.") +
geom_line(aes(colour=(ORF), group=(Colony)), linetype=2) +
scale_colour_manual(values=c("black", "gray"),
guide = guide_legend(title.position = "top", nrow = 2),
name = "Pocillopora type",
labels=c("Type 3", "Type 1")) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.01,0.01)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
Traject_ProportionORF <- Traject_ProportionORF_a +
geom_jitter((aes(colour=factor(ORF))),
height = 0, width = 0.3, color="black", size=1, shape=21) +
scale_fill_manual(values=c("white", "gray"))
#Traject_ProportionORF
Figure_2_ORF<-grid.arrange(Histo2014_DpropORF, Traject_ProportionORF, Histo2016_DpropORF,
nrow=1, widths=c(3/9, 3/9, 3/9))
# ggsave(file="Outputs/Figure_2_ORF.svg",plot=Figure_2_ORF, width=7, height=3.8, dpi=300)
data.Uva$Mor_ITS2<-paste(data.Uva$DominantITS2, data.Uva$Mortality, sep = "_")
data.Uva$Mor_ITS2<-factor(data.Uva$Mor_ITS2, levels = c(
"C1d-dominated_No", "C1d-dominated_Partial", "C1d-dominated_Almost total",
"C1bc-dominated_No", "C1bc-dominated_Partial",
"D1-dominated_No", "D1-dominated_Partial", "D1-dominated_Almost total"))
summary(data.Uva$Mor_ITS2)
## C1d-dominated_No C1d-dominated_Partial
## 11 3
## C1d-dominated_Almost total C1bc-dominated_No
## 4 17
## C1bc-dominated_Partial D1-dominated_No
## 3 43
## D1-dominated_Partial D1-dominated_Almost total
## 3 3
Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(Mor_ITS2))) +
ggthe_bw + coord_flip()+ ggtitle("a.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("skyblue1",
"darkorchid1",
"red1"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (August 2014)") +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~((D/H)/(S/H))),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2014_DpropITS2
Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(Mor_ITS2))) +
ggthe_bw + coord_flip()+ ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("skyblue1", "blue2", "blue4",
"darkorchid3",
"red1", "red3", "red4"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (April 2016)") +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Histo2016_DpropITS2
HistoAll_DpropITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(Mor_ITS2))) +
ggthe_bw + coord_flip()+ ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("skyblue1", "blue2", "blue4",
"darkorchid1","darkorchid3",
"red1", "red3", "red4"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (April 2016)") +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA)) + facet_grid(ORF~Year)
#HistoAll_DpropITS2
Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp, group=(Colony)))+
ggthe_bw + ggtitle("b.") +
geom_line(aes(colour=ITS2_14),linetype=2) +
geom_jitter(aes(colour=Mor_ITS2, shape=Mortality),
height = 0.00, width = 0.25, alpha=0.9) +
scale_colour_manual(values=c(
"darkorchid1", "darkorchid1","darkorchid3",
"skyblue1","blue4", "skyblue1", "blue2",
"red1", "red4", "red1", "red3" ),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (April 2016)") +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.01,0.01)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))
#Traject_ProportionITS2
Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(3/9, 3/9, 3/9))
# ggsave(file="Outputs/Figure_ITS2_Mortality.svg", plot=Figure_2_ITS2, width=7, height=4)
data.Uva$Mor_ITS2<-paste(data.Uva$DominantITS2, data.Uva$Mortality, sep = "_")
data.Uva$Mor_ITS2<-factor(data.Uva$Mor_ITS2, levels = c(
"C1d-dominated_Almost total", "C1d-dominated_Partial","C1d-dominated_No",
"C1bc-dominated_Partial","C1bc-dominated_No",
"D1-dominated_Almost total", "D1-dominated_Partial", "D1-dominated_No" ))
summary(data.Uva$Mor_ITS2)
## C1d-dominated_Almost total C1d-dominated_Partial
## 4 3
## C1d-dominated_No C1bc-dominated_Partial
## 11 3
## C1bc-dominated_No D1-dominated_Almost total
## 17 3
## D1-dominated_Partial D1-dominated_No
## 3 43
Histo2014_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2014), ], aes (D.Prp , fill=(Mor_ITS2))) +
ggthe_bw + ggtitle("a.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("black",
"black",
"black"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (August 2014)") +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~((D/H)/(S/H))),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2014)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="none",
legend.title.align=0.5,
strip.background =element_rect(fill=NA)) + facet_grid(ORF~.)
#Histo2014_DpropITS2
Histo2016_DpropITS2 <- ggplot(data.Uva[which(data.Uva$Year.2==2016), ], aes (D.Prp , fill=(Mor_ITS2))) +
ggthe_bw + ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("white", "gray", "black",
"gray",
"white", "gray", "black"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (April 2016)") +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="none",
legend.title.align=0.5,
strip.background =element_rect(fill=NA)) + facet_grid(ORF~.)
#Histo2016_DpropITS2
HistoAll_DpropITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(Mor_ITS2))) +
ggthe_bw + ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0) +
scale_fill_manual(values=c("skyblue1", "blue2", "blue4",
"darkorchid1","darkorchid3",
"red1", "red3", "red4"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (April 2016)") +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies \n (Aug 2016)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA)) + facet_grid(ORF~Year)
#HistoAll_DpropITS2
HistoAll_DpropITS2 <- ggplot(data.Uva, aes (D.Prp , fill=(Mor_ITS2))) +
ggthe_bw + ggtitle("c.")+
geom_histogram(binwidth = 0.05, boundary=0, color = "black") +
scale_fill_manual(values=c("white", "gray", "black",
"gray","black",
"white", "gray", "black"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 community (April 2016)") +
scale_x_continuous(name="",
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies"),
limits = c(0,14),
breaks=seq(0,14,by=2)) +
theme(legend.position="bottom",
legend.title.align=0.5,
strip.background =element_rect(fill=NA)) +
facet_grid(ORF~Year) + theme(legend.position = "none")
#HistoAll_DpropITS2
#ggsave(file="Outputs/Panels.svg", plot=HistoAll_DpropITS2, width=6, height=2.8)
Traject_ProportionITS2<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp, group=(Colony)))+
ggthe_bw + ggtitle("b.") +
geom_line(linetype=2) +
geom_jitter(aes(colour=Mortality, shape=ORF), size=1.5,
height = 0.00, width = 0.25, alpha=1) +
scale_colour_manual(values=c("black", "gray","white")) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.02,0.02)) +
theme(legend.position="none",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))# + facet_grid(ORF~.)
#Traject_ProportionITS2
#ggsave(file="Outputs/Trajectory.svg", plot=Traject_ProportionITS2, width=6.5, height=3)
Figure_2_ITS2<-grid.arrange(Histo2014_DpropITS2, Traject_ProportionITS2,Histo2016_DpropITS2,
nrow=1, widths=c(2.5/9, 4/9, 2.5/9))
# ggsave(file="Outputs/Figure_ITS2_Mortality2.svg", plot=Figure_2_ITS2, width=7, height=4)
DpropORF <- ggplot(data.Uva, aes (D.Prp , fill=(Mortality))) +
ggthe_bw + ggtitle("a.")+
geom_histogram(binwidth = 0.1, boundary=0, color="black") +
scale_fill_manual(values=c("black",
"gray",
"white"),
guide = guide_legend(title.position = "top", nrow = 1),
name = "Colony condition") +
scale_x_continuous(name= expression(
italic(Durusdinium)~proportion~((D/H)/(S/H))),
breaks=(seq(0,1,by=0.1)),
expand=c(0.01,0.01)) +
scale_y_continuous((name= "Number of colonies
\n (Aug 2014)
\n (Aug 2014) \n (Aug 2014)"),
limits = c(0,17),
breaks=seq(0,15,by=5)) +
theme(legend.position="top",
legend.title.align=0.5,
strip.background =element_rect(fill=NA)) + facet_grid(ORF~Year)
#DpropORF
Traject_ORF<- ggplot(data.Uva, aes (x=Year.2, y=D.Prp, group=(Colony)))+
ggthe_bw + ggtitle("b.") +
geom_line(linetype=2) +
geom_jitter(aes(colour=Mortality, shape=ORF), size=2,
height = 0.00, width = 0.25, alpha=1) +
scale_colour_manual(values=c("black", "gray","white")) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016")) +
scale_y_continuous((name= ""),
breaks=seq(0,1,by=0.1),
expand=c(0.02,0.02)) +
theme(legend.position="none",
legend.title.align=0.5,
strip.background =element_rect(fill=NA))# + facet_grid(ORF~.)
#Traject_ORF
Figure_3bc_ORF<-grid.arrange(Traject_ORF, DpropORF,
nrow=2, heights=c(1.2/3, 1.8/3))
#ggsave(file="Outputs/Fig3_v8.svg", plot=Figure_3bc_ORF, width=6.5, height=7)
# Summary dominant genera
Community_prescence_absence<-data.Uva %>%
group_by(ORF, Community, Year.2) %>%
dplyr::summarise(total.count=n())
Community_prescence_absence
Community_dominance<-data.Uva %>%
group_by(ORF, DominantCommunity, Year.2) %>%
dplyr::summarise(total.count=n())
Community_dominance
Community_dominance<-data.Uva %>%
group_by(ORF, Year.2, DominantCommunity,
Community) %>%
dplyr::summarise(total.count=n())
Community_dominance
# SH_mean<-data.Uva %>%
# group_by(ORF, Year.2, DominantCommunity) %>%
# dplyr::summarise(mean=mean(tot.SH))
# SH_mean
#
# SH_SE<-data.Uva %>%
# group_by(ORF, Year.2, DominantCommunity) %>%
# dplyr::summarise(sd=sd(tot.SH))
# SH_SE
#SH_2014CoI<-lm(logTot.SH ~ Year * ORF * DominantCommunity, data=data.Uva)
SH_2014CoI<-lmerTest::lmer(logTot.SH ~ Year * ORF * DominantCommunity + (1|Colony), data=data.Uva)
summary(SH_2014CoI)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logTot.SH ~ Year * ORF * DominantCommunity + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 75.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3367 -0.5185 -0.0403 0.3713 3.4413
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.01992 0.1411
## Residual 0.10475 0.3236
## Number of obs: 87, groups: Colony, 29
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -1.68521 0.11691 77.39879
## Year15 0.22667 0.15257 51.72506
## Year16 -0.05117 0.26368 68.09623
## ORFType_3 0.02507 0.18560 75.84196
## DominantCommunityDurusdinium-dominated 0.14102 0.14919 77.99999
## Year15:ORFType_3 -0.79408 0.24123 51.72506
## Year16:ORFType_3 -0.72686 0.32318 62.79186
## Year15:DominantCommunityDurusdinium-dominated -0.06059 0.19555 51.72506
## Year16:DominantCommunityDurusdinium-dominated 0.01196 0.28860 68.27239
## t value Pr(>|t|)
## (Intercept) -14.415 <2e-16 ***
## Year15 1.486 0.1434
## Year16 -0.194 0.8467
## ORFType_3 0.135 0.8929
## DominantCommunityDurusdinium-dominated 0.945 0.3475
## Year15:ORFType_3 -3.292 0.0018 **
## Year16:ORFType_3 -2.249 0.0280 *
## Year15:DominantCommunityDurusdinium-dominated -0.310 0.7579
## Year16:DominantCommunityDurusdinium-dominated 0.041 0.9671
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 ORFT_3 DmnCD- Y15:OR Y16:OR Y15:DC
## Year15 -0.653
## Year16 -0.373 0.289
## ORFType_3 -0.630 0.411 0.235
## DmnntCmmnD- -0.777 0.511 0.290 0.489
## Yr15:ORFT_3 0.413 -0.632 -0.183 -0.650 -0.323
## Yr16:ORFT_3 0.305 -0.236 -0.816 -0.483 -0.237 0.373
## Yr15:DmnCD- 0.509 -0.780 -0.226 -0.321 -0.655 0.493 0.184
## Yr16:DmnCD- 0.359 -0.264 -0.920 -0.226 -0.463 0.167 0.750 0.339
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients
anova(SH_2014CoI)
ranova(SH_2014CoI)
step(SH_2014CoI, reduce.random=FALSE)
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -37.966 97.933
## (1 | Colony) 0 10 -38.884 97.768 1.8355 1 0.1755
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Year:ORF:DominantCommunity 0
##
## Model found:
## logTot.SH ~ Year * ORF * DominantCommunity + (1 | Colony)
SH_2014CoI<-lm(logTot.SH ~ Year * ORF * DominantCommunity, data=data.Uva)
SH_I.emmc<-emmeans(SH_2014CoI, ~Year*ORF*DominantCommunity)
SH_I_groups<-multcomp::cld(SH_I.emmc)
SH_I_groups<-as.data.frame(SH_I_groups[complete.cases(SH_I_groups),])
SH_I_groups<-SH_I_groups[order(
SH_I_groups$Year, SH_I_groups$ORF),]
SH_I_groups
#write.csv(SH_I_groups,"Outputs/SH_ORF_Dominant_groups2.csv")
# plot(emmeans(SH_2014CoI, ~Year|~DominantCommunity), comparisons = TRUE) + coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_wrap(~DominantITS2)+ ggthe_bw
SH_2014ComI <- ggplot(data.Uva, aes (Year.2, logTot.SH, colour=factor(DominantITS2))) + ggthe_bw + facet_wrap(~ORF) + ggtitle("a.")+
#geom_point(alpha=0.5)+
#geom_line(aes(group=(Colony)))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2,
position=position_dodge(width=0.5))+
stat_summary(fun.data = "mean_cl_boot",geom = "point",
position=position_dodge(width=0.5))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.5)) +
scale_colour_manual(values=c("blue3","blue3" , "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant symbiont",
labels=c("Cladocopium", "Durusdinium"))+
theme(legend.position="mone",
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-3, -0.1),
name="\n Total symbiont to host cell ratio (S/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" ", "", " "))
SH_2014ComI
#ggsave(file="Figure_3.svg", plot=SH_ITS2_Comunities, width=6, height=4, dpi=300)
CH_Com2ITS2 <- ggplot(data.Uva, aes (Year.2, logC.SH, colour=factor(D14Dominant))) + ggthe_bw+ ggtitle("c.") + facet_grid(~ORF) + #geom_point(alpha=0.5)+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue3", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant symbiont",
labels=c("Cladocopium", "Durusdinium"))+
theme(legend.position="mone",
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-6, -0.1),
name="\n Cladocopium to host cell ratio (C/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
labels = c(" ", " ", " "))
CH_Com2ITS2
CH_I<-lmerTest::lmer(logC.SH ~ Year* DominantITS2 + (1|Colony), data=data.Uva)
summary(CH_I)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logC.SH ~ Year * DominantITS2 + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 108.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.89915 -0.40700 -0.07625 0.38592 1.92091
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.06891 0.2625
## Residual 0.38254 0.6185
## Number of obs: 55, groups: Colony, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.68510 0.22354 44.93743 -7.538 1.65e-09
## Year15 0.22932 0.29156 28.78344 0.787 0.437988
## Year16 -0.12883 0.50336 39.26788 -0.256 0.799334
## DominantITS2C1d-dominated 0.02496 0.35386 44.57426 0.071 0.944077
## DominantITS2D1-dominated -1.27282 0.35244 45.62565 -3.611 0.000755
## Year15:DominantITS2C1d-dominated -0.79673 0.46100 28.78344 -1.728 0.094658
## Year16:DominantITS2C1d-dominated -0.64921 0.61716 35.78277 -1.052 0.299886
## Year15:DominantITS2D1-dominated -0.44588 0.49919 31.70568 -0.893 0.378478
## Year16:DominantITS2D1-dominated 0.31093 0.62644 42.90497 0.496 0.622184
##
## (Intercept) ***
## Year15
## Year16
## DominantITS2C1d-dominated
## DominantITS2D1-dominated ***
## Year15:DominantITS2C1d-dominated .
## Year16:DominantITS2C1d-dominated
## Year15:DominantITS2D1-dominated
## Year16:DominantITS2D1-dominated
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 DITS2C DITS2D Y15:DITS2C Y16:DITS2C Y15:DITS2D
## Year15 -0.652
## Year16 -0.376 0.290
## DmnnITS2C1- -0.632 0.412 0.238
## DmnnITS2D1- -0.631 0.414 0.238 0.398
## Y15:DITS2C1 0.412 -0.632 -0.183 -0.651 -0.262
## Y16:DITS2C1 0.307 -0.236 -0.816 -0.486 -0.194 0.373
## Y15:DITS2D1 0.382 -0.584 -0.170 -0.241 -0.606 0.369 0.138
## Y16:DITS2D1 0.339 -0.233 -0.816 -0.214 -0.538 0.147 0.665 0.345
anova (CH_I)
#step(CH_I)
CH_I.emmc<-emmeans(CH_I, ~Year*DominantITS2)
CH_I_groups<-cld(CH_I.emmc)
CH_I_groups<-as.data.frame(CH_I_groups[complete.cases(CH_I_groups),])
CH_I_groups<-CH_I_groups[order(CH_I_groups$DominantITS2,CH_I_groups$Year),]
CH_I_groups
#write.csv(CH_I_groups,"CH_I_groups.csv")
DH_Com2ITS2 <- ggplot(data.Uva, aes (Year.2, logD.SH, colour=factor(D14Dominant))) + # geom_point(alpha=0.5)+
ggthe_bw+ ggtitle("d.")+ facet_grid(~ORF) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, position=position_dodge(width=0.3))+
stat_summary(fun.data = "mean_cl_boot",geom = "point", position=position_dodge(width=0.3))+
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=0.3)) +
scale_colour_manual(values=c("blue", "red3"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant symbionts",
labels=c("Cladocopium-dominated", "Durusdinium-dominated"))+
theme(legend.position=c(0.75,0.5),
legend.box.background = element_rect(),
strip.background =element_rect(fill="white"))+
scale_y_continuous(limits = c(-6, -0.1),
name="\n Durusdinium to host cell ratio (D/H)",
labels = (math_format(10^.x)),
expand = c(0.0, 0)) +
scale_x_continuous(name="",
breaks = c(2014, 2015, 2016),
# labels = c(" Aug \n 2014", " Aug \n 2015", " Apr \n 2016"))
labels = c(" ", " ", " "))
DH_Com2ITS2
DH_LM_ITS2<-lmer(logD.SH ~ Year* DominantITS2 + (1|Colony), data=data.Uva)
summary(DH_LM_ITS2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logD.SH ~ Year * DominantITS2 + (1 | Colony)
## Data: data.Uva
##
## REML criterion at convergence: 62.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1646 -0.5042 0.0345 0.3558 3.7511
##
## Random effects:
## Groups Name Variance Std.Dev.
## Colony (Intercept) 0.000 0.0000
## Residual 0.158 0.3975
## Number of obs: 58, groups: Colony, 22
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -4.86573 0.16229 52.00000 -29.981 < 2e-16
## Year15 0.07049 0.32459 52.00000 0.217 0.829
## Year16 2.95052 0.42939 52.00000 6.871 7.93e-09
## DominantITS2D1-dominated 3.26614 0.19398 52.00000 16.838 < 2e-16
## Year15:DominantITS2D1-dominated 0.14077 0.35768 52.00000 0.394 0.696
## Year16:DominantITS2D1-dominated -2.96086 0.45076 52.00000 -6.569 2.41e-08
##
## (Intercept) ***
## Year15
## Year16 ***
## DominantITS2D1-dominated ***
## Year15:DominantITS2D1-dominated
## Year16:DominantITS2D1-dominated ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year15 Year16 DITS2D Y15:DI
## Year15 -0.500
## Year16 -0.378 0.189
## DmnnITS2D1- -0.837 0.418 0.316
## Y15:DITS2D1 0.454 -0.907 -0.171 -0.542
## Y16:DITS2D1 0.360 -0.180 -0.953 -0.430 0.233
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#step(DH_LM_ITS2)
DH_I.emmc<-emmeans(DH_LM_ITS2, ~Year*DominantITS2)
DH_groups_I<-cld(DH_I.emmc)
DH_groups_I<-as.data.frame(DH_groups_I[complete.cases(DH_groups_I),])
DH_groups_I<-DH_groups_I[order(DH_groups_I$DominantITS2,DH_groups_I$Year),]
DH_groups_I
#write.csv(DH_groups_I,"Output/DH_I_groups.csv")
Figure_4<-grid.arrange(SH_2014ComI, CH_Com2ITS2, DH_Com2ITS2,nrow=3)
Figure_4_InitailCommunity<-grid.arrange(CH_Com2ITS2, DH_Com2ITS2,nrow=2)
# ggsave(file="Outputs/Figure_3c.svg", plot=Figure_3, width=7, height=11)
# ggsave(file="Outputs/Figure_3d.png", plot=Figure_3_InitailCommunity, width=7, height=7)
Select data
log_Log.data<-data.Uva[,c("Colony","Year.2", "C.SH", "D.SH",
"D14Community", "Community","ITS2_14",
"ORF", "DominantITS2", "Mortality", "Mor_ITS2")]
log_Log.data<-log_Log.data %>%
filter(Year.2 != "2015")
log_Log.data$Yr<-as.factor(log_Log.data$Year.2)
log_Log.data$Yr_M<-as.factor(
paste(log_Log.data$Yr,
log_Log.data$Mortality,
sep = "_"))
log_Log.data$logC<-log10(log_Log.data$C.SH)
log_Log.data$logD<-log10(log_Log.data$D.SH)
log_Log.data$logC[which(log_Log.data$C.SH==0)] <- -6
log_Log.data$logD[which(log_Log.data$D.SH==0)] <- -6
Community_Changes1 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=ITS2_14, shape=ORF)) +
geom_abline(intercept=0, slope=1, linetype=1, color="black")+
geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
geom_abline(intercept=1, slope=1, linetype=2, color="black")+
geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Dominant ITS2 types (2014)",
labels=c("C1d", "C1b-c", "D1"))+
ggthe_bw + theme(legend.position="bottom",
legend.box = "horizontal")+
scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
geom_point(alpha=0.5, size=4) +
geom_path(size = 1,arrow =
arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
"10,0000D : 1C", "1,000D : 1C", "100D : 1C","10D : 1C", "C : D"), angle = 45)+
annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
"1D : 10,000C", "1D : 1,000C", "1D : 100C", "1D : 10C"), angle = 45)
Community_Changes1
#ggsave(file="FigAndrewsFavoriteGraph_ITS2.png", plot=Community_Changes1, width=5.5, height=7.1)
Community_Changes2 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=ITS2_14)) +
geom_abline(intercept=0, slope=1, linetype=1, color="black")+
geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
geom_abline(intercept=1, slope=1, linetype=2, color="black")+
geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
scale_colour_manual(values=c("blue","purple", "red"),
guide = guide_legend(title.position = "top", nrow = 3),
name = "Initial ITS2",
labels=c("Cd1", "C1b-c", "D1"))+
#scale_colour_manual(values=c("skyblue1", "blue2", "blue4",
# "darkorchid1","darkorchid3",
# "red1", "red3", "red4"))+
ggthe_bw +
theme(legend.position="bottom",
legend.box = "vertical")+
#theme(legend.position="none")+
scale_x_continuous(limits = c(-6,-0.1),
"Cladocopium abundance (log10 C/H))")+
scale_y_continuous(limits = c(-6,-0.1),
"Durusdinium abundance (log10 D/H))") +
geom_point(aes(shape=Mortality),
alpha=0.5, size=3) +
geom_path(size = 0.5 , arrow =
arrow(type = "open", angle = 30,
length = unit(0.1, "inches")))+
annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
"10,0000D : 1C", "1,000D : 1C", "100D : 1C","10D : 1C", "C : D"), angle = 45)+
annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
"1D : 10,000C", "1D : 1,000C", "1D : 100C", "1D : 10C"), angle = 45)
Community_Changes2
#ggsave(file="LogLog_labels.svg", plot=Community_Changes2, width=5.5, height=6)
Community_Changes3 <- ggplot (log_Log.data, aes(logC, logD, gropup=Colony, colour=Colony)) +
geom_abline(intercept=0, slope=1, linetype=1, color="black")+
geom_abline(intercept=-1, slope=1, linetype=2, color="black")+
geom_abline(intercept=-2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=-3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=-5, slope=1, linetype=2, color="gray")+
geom_abline(intercept=1, slope=1, linetype=2, color="black")+
geom_abline(intercept=2, slope=1, linetype=2, color="darkgray")+
geom_abline(intercept=3, slope=1, linetype=2, color="gray")+
geom_abline(intercept=4, slope=1, linetype=2, color="gray")+
geom_abline(intercept=5, slope=1, linetype=2, color="gray")+
#scale_colour_manual(values=c("blue","purple", "red"),
# guide = guide_legend(title.position = "top", nrow = 3),
# name = "Detected symbionts (2014)",
# labels=c("Cladocopium", "Mixed", "Durusdinium"))+
ggthe_bw + theme(legend.position="none")+
scale_x_continuous(limits = c(-6,-0.1), "log10 (C/H cell ratio)")+
scale_y_continuous(limits = c(-6,-0.1), "log10 (D/H cell ratio)") +
geom_point(aes(shape=Yr), alpha=0.5, size=4) +
geom_path(size = 1,arrow =
arrow(type = "open", angle = 30, length = unit(0.1, "inches")))+
annotate("text", x = c(seq(-5.1, -1.1, by=1)), y = -0.9, label = c(
"10,0000D : 1C", "1,000D : 1C", "100D : 1C","10D : 1C", "C : D"), angle = 45)+
annotate("text", x = -0.9, y = c(seq(-4.7, -1.7, by=1)), label = c(
"1D : 10,000C", "1D : 1,000C", "1D : 100C", "1D : 10C"), angle = 45)
Community_Changes3
# ggsave(file="LogLog_Colony.svg", plot=Community_Changes3, width=6, height=6)
# Only type 1 colonies (n=23)
Input =("
Year C D
2014 9 14
2016 2 18
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.03936
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.9354645 61.0268633
## sample estimates:
## odds ratio
## 5.560207
# Type 1 and type 3 colonies
Input =("
Year C D
2014 15 14
2016 4 18
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.01988
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.145585 23.762126
## sample estimates:
## odds ratio
## 4.670456
# Spp changes
Input =("
Year I II
2014 23 6
2016 20 2
")
Matriz <- as.matrix(read.table(textConnection(Input),
header=TRUE,
row.names=1))
fisher.test(Matriz,
alternative="two.sided")
##
## Fisher's Exact Test for Count Data
##
## data: Matriz
## p-value = 0.4399
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.03477518 2.50826671
## sample estimates:
## odds ratio
## 0.3901317